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Minimal Boltzmann kinetic models, such as lattice Boltzmann, are often used as an 
alternative to the discretization of the Navier-Stokes equations for hydrodynamic simu- 
lations. Recently, it was argued that modeling sub-grid scale phenomena at the kinetic 
level might provide an efficient tool for large scale simulations. Indeed, a particular 
variant of this approach, known as the entropic lattice Boltzmann method (ELBM), has 
shown that an efficient coarse-grained simulation of decaying turbulence is possible us- 
ing these approaches. The present work investigates the efficiency of the entropic lattice 
Boltzmann in describing flows of engineering interest. In order to do so, we have chosen 
the flow past a square cylinder, which is a simple model of such flows. We will show 
that ELBM can quantitatively capture the variation of vortex shedding frequency as a 
function of Reynolds number in the low as well as the high Reynolds number regime, 
without any need for explicit sub-grid scale modeling. This extends the previous studies 
for this set-up, where experimental behavior ranging from Re ~ O(10) to Re < 1000 
were predicted by a single simulation algorithm I I 1 ^ 1 1 ^ 1 

Keywords: Lattice Boltzmann; H theorem; Turbulence Model; Flow Past Square Cylin- 
der 

1. Introduction 

Direct numerical simulations of low to moderate Reynolds number turbulent flows 
in simple geometries, such as channel or infinite domain, have substantially im- 
proved the understanding of turbulent flows. However, for fully developed turbulent 
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flows of engineering interest, such simulations are far beyond the processing capa- 
bilities of even the most powerful computers. Due to the huge difference between 
the smallest length scale, Ik, (known as Kolmogorov length) and the largest length 
scale (typically on the scale of system-size), all length scales cannot be taken into 
account in any simulations. On the other hand, the convective non-linearity in the 
Navier-Stokes equations ensures that hydrodynamic equation at all length scales 
are coupled to each other. Further complications arise due to the non-linearity in- 
troduced through the pressure term and the boundary conditions. Thus, writing 
model equations for large scale turbulence simulations, which can predict the ef- 
fects of unresolved (small) scales of motion on the resolved (large) ones constitutes 
a major part of turbulence research. 

Many models for turbulence rely on the idea of scale-separation, which assumes 
that there exist a cut-off length scale, l c , (l c ^> Ik) such that small-scale eddies, 
I < l c , are in local equilibrium with large scale eddies (/ > l c ). For example, all eddy- 
viscosity models rely on the analogy between turbulent and molecular transport. 
Within this analogy, the small-scale eddies play the role of the molecules, while the 
Kolmogorov length-scale Ik, (in actual computations replaced by the grid spacing 
S) plays the role of the mean-free path. This analogy with kinetic theory has the 
significant advantage of leaving the form of the Navier-Stokes equation for the large 
eddies intact, with an effective viscosity v e replacing the molecular viscosity v. 

However, the scale-separation argument discussed above simply does not hold 
for turbulence. Thus, these models often fail to reproduce strongly off-equilibrium 
turbulent effects. In particular, the scenario described by the eddy viscosity model 
clearly indicates that the dynamics of the smallest resolved scales (I — > l c ) faces a 
situation similar to finite-Knudsen flows at Kn ~ 1. This lack of scale separation 
is the fundamental reason for the failure to develop a consistent theory of fully 
developed turbulence. For the classical Boltzmann equation, the proble m of Kn ~ 1 
hydrodynamics was solved by the method of invariant manifolds 1 ^ 1 718 1 

It has been recently pointed out that the kinetic representation of hydrody- 
namics provides a natural generalization of the notion of eddy viscosity to such 
non-equilibrium high-Kn t regimes | 9 | 10 | ll|l ll (Kn t — Ik /I, for an eddy of size I). 
The key point is that solutions of the kinetic equation apply at all orders of the 
Knudsen number, so that any kinetic model ensuring correct hydrodynamic behav- 
ior would handle the dynamics of small eddies at Kn t ~ 1 in a way which goes 
beyond the eddy viscosity representation. 

Further, important computational features of the kinetic models are the linear- 
ity of the streaming operator and the locality (in configuration space) of the the 
non-linearity introduced by the collision term. In order to efficiently use these ad- 
vantages of kinetic models, while getting rid of enormously large numbers of degree 
of freedoms associated with the (true) Boltzmann equation, a drastically simplified 
versions of the Boltzmann equation, known as the lattice Boltzmann method^2^31 
(hereafter LBM), had been developed in the last decade. Indeed, the m ethod has 
been applied to a large variety of fluid flows, including turbulent ones | 12 | 15 | 16 | 
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However, owing to the lack of a H-theorem 021 the standard lattice Boltzmann 
method often exhibits disruptive non-linear instabilities associated with the dy- 
namics of near-grid scales with I ~ S. In order to improve on this shortcoming 
of the method, while preserving its computational efficiency, a modification of the 
algorithm, known as "entropic lattice Boltzmann method" , was developed recently 
| 18 I ^U I ^1 | ^3J25J The compliance with the H -theorem ensures the non-linear sta- 
bility for the method. The basic idea is to exploit the requirement of increase in 
the entropy ( "grid entropy" ) as a criterion to decouple the resolved scales from the 
unresolved ones. The compliance with the H-theorem ensures that the termination 
of the cascade picture on the grid scale does not lead to non-physical disruptive 
non-linear instabilities. 

In fact, the model with its adaptive but local regulation of the relaxation time, 
can be regarded as a turbulence model inspired by genuinely kinetic requirements. 
The efficiency of suc h tu rbulence modeling was recently tested in the context of the 
decaying turbulence 1^41. In this work, we test the efficiency of the method for the 
flow past a square cylinder, a prototype of flows of engineering interest. This set-up 
is used often to test the efficiency of different Navier-Stokes (or its model) solvers 
(for example see Refs. 1-5). Even though the set-up is simple enough, experimen- 
tal results predict a n on-trivial dependence of the vortex-shedding frequency on the 
Reynolds number l^fil. The vortex-shedding frequency is characterized by a dimen- 
sionless number known as Strouhal number, which is defined as S = nLc^/Uo, 
where n is the frequency of the vortex shedding, L c h a r is the length of the cylinder, 
and Uq is the characteristic velocity, taken as the velocity at the inlet. The Strouhal 
number undergoes a sharp transition around Re ~ 70 and attends a more or less 
constant value of 0.13 at high Reynolds numbers EEEZI Although the transition at 
Re ~ 70 is captured by different schemes with varying degree of accuracy, so far 
to the best of our knowledge using direct Numerical simulation of Navier-Stokes 
equation, the experimental curve is reproduced up to Re = 1000. Various models of 
turbulence have also be en tested on the flow past a square cylinder at Re = 21400 
pg | ^ | 3U I 31 | 32 | 33 | ;j4 | 3g ) We wi u sll0w tlmt ELBM can quantitatively capture the 

variation of vortex shedding frequency as a function of Reynolds number without 
any need for explicit sub-grid scale modeling. 

The work is organized as follows: A brief description of the LBM and its entropic 
version is provided in Sec. II. In Sec. Ill, implementation details for the flow past a 
square cylinder will be provided. In Sec. IV, a comparison of the present work with 
the experimental results and other numerical results will be presented. Further, in 
Sec. V, the outlook for the present work will be discussed. 



2. The entropic Lattice Boltzmann method 



In the lattice Boltzmann method, the discrete kinetic equation for the one particle 
distribution function, fi = f(x,Ci,t), which denotes the probability of finding a 
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particle with velocity Cj at position x and time t, is written as: 

dtfi + C< ■ dxfi = -T- 1 (/« - Jf) , (1) 

where, the right-hand-side of this equation represents collisional relaxation to the 
local equilibrium, /° q , on a time scale r. The irreversible relaxation to this local 
equilibrium provides viscous behavior, with a kinematic viscosity of the order of 
v ~ c 2 s t, where c s is the speed of sound. A general formulation of the kinetic 
models is presented in the work of Gorban and Karlin 1 ^ 1 . 

The set of discrete velocities is chosen in such a way as to ensure sufficient sym- 
metry to recover the conservation of mass, momentum and momentum-flux. The 
minimal set of discrete velocities needed to reconstruct the Navier-Stokes equations 
(as the large-scale limit) are related to zeroes of third-order Hermite polynomials 
in Cj. For example in one-dimension (spatial dimension D = 1), the three discrete 
velocities are 



C = c {-1,0,1}, (2) 

where c is some arbitrary constant usually taken as c = Sx/St, with St and Sx being 
step-size in the spatial and the temporal discretization, respectively. In higher 
dimensions, the discrete velocities are tensor products of the discrete velocities in 
one dimension. 

The discrete-velocity local equilibrium is the minimizer of the the discrete H 
function: 

# W ,c l} =£/an(A) (3 ) 

i— 1 



where, Wj are the weights associated with each discrete velocity (for the one- 
dimensional case) 

»'=aui (4) 



,6' 3' 6, 

under the constraint of local conservation laws: 

3 d 

^/rU, C,;} = {p, P U}. (5) 

i=l 

In higher dimensions, the weights are constructed by multiplying weights associated 
with each direction. 

The explicit expression for f° q is^Q 



with j being the index for spatial directions. 
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In order to discretise the model kinetic equation in a way consistent with the 
second law of thermodynamics (_ff-theorem) , the notion of the bare collision A, de- 
fined as the collision term stripped of its relaxation parameters, is introduced. For 
the case of the single relaxation time collision model considered here (right hand 
side of Eq. [JJ, the bare collision is given as A = f cq — f, where f denotes the 
S^-dimensional population vector. The time stepping in this method is performed 
through an over-relaxation collisional process and linear convection through a se- 
quence of steps in which the H function is bound not to decrease. The monotonicity 
constraint on the H function is imposed through the following geometric procedure: 
In the first step, populations are changed in the direction of the bare collision in 
such a way that the H function remains constant. In the second step, dissipation 
is introduced and the magnitude of the H function decreases. Thus, 



/ i (x,^) = /i(x-C i ^,0) + a/3 



f^(x.-CiSt,Q)-f(it.-C i 5t,0) 



(7) 



were j3 is the discrete form of the relaxation frequency related to r as 



and the parameter a is defined by the condition: 

H (f ) = H (f + oA) . (9) 

Close t o th e local equilibrium a is equal to 2 (Implementation details can be found 
m Ref. H3). Equation!! is the essence of the discrete time iJ-theorem for the discrete 
kinetic equation The notion of the discrete time H theorem was first introduced 
by Karlin et al^Sl It is the ext ension of the earlier general study of the initial layer 
problem in dissipative kinetics ESEDI The local adjustments of the relaxation time 
(via the parameter a), as dictated by compliance with the H theorem, guarantee 
positivity of the distribution function also for the case of discrete steps, thereby 
ensuring non-linear stability of the numerical scheme. 

In this method, for the solid-fluid boundary condition a discret izat ion scheme of 
diffusive boundary condition for the Boltzmann equation is used^2 The essence 
of the diffusive boundary condition is that particles loose their memory of the in- 
coming direction after reaching the wall. Once a particle reaches the wall, it gets 
redistributed in a way consistent with the mass-balance and normal-flux condi- 
tions. Further, the boundary condition must also satisfy the condition of detailed 
balance: if the incoming populations are at equilibrium (corresponding to the wall- 
velocity), the outgoing populations are also at equilibrium (corresponding to the 
wall- velocity) . As an example, let us consider the case when the wall normal, e, 
(pointing towards the fluid) is in the positive y direction (see Fig. For this 
particular case, the boundary update rules for incoming and grazing populations 
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Figure 1: A schematic diagram for situation near a fiat wall, when the wall normal, 
e, (pointing towards the fluid) is in the positive y direction. Here, the case of 
two-dimensional D2Q9 lattice is considered. 



on a two-dimensional lattice are: 

fo(x,y,t + St) = fo(x,y,t), 
fi(x,y,t + St) = fo(x - cSt,y,t), 
f 2 {x,y,t + St) = fo(x + cSt,y,t), 

U(x, y,t + 6t) = ± \St (x, y + cSt, t) + ft (x, y, *)] , (10) 

f 7 (x,y,t+ St) = ~ [/ 7 * (x + cSt, y + cSt, t) + f* 7 (x, y, *)] , 

fs(x,y,t + St) = - [/|(x - cSt,y + cSt,t) + fg(x,y,t)] , 

where /* denotes after collision populations, and update rules for outgoing popu- 
lations are: 



wall ) 



_ „ , / 4 (x, t + St) + / 7 (x, t + St) + / 8 (x, t + St) 

f2 q (P: U wall) + /^(p, U wan ) + f^{p, U wal i) 

f( +mx+\ M \ /4(x ' 1 ± St) ± /?(x ' t + 8t) + fs(x,t + St) 

f 5 (x,t + 5t)=f 5 ( P , iw) ^ — } - + /r(p; uwau) (11) 

h (Pi U wall) + fsiP, U W all) + / 6 (P, "wall) 



3. Implementation Details for Flow Past a Square Cylinder 

We study the flow past a square cylinder in the two-dimensional set-up shown 
schematically in Fig. [21 The square cylinder is placed at a location L\ = 10 L c har 
downstream form the inlet and symmetric in the other direction. For high Reynolds 
number simulations (above Re > 1000), a computational grid of width W = 25L c i iar 
and length L = 45 L c har was used, while for moderate Reynolds number simulations 
(above Re < 1000), a computational grid of width W — 25L c i iar and length L = 
30 i c har was used. The width of the domain results in a blockage ratio (£ c har/ W) 
of less than 5. In all simulations the inlet velocity J7 is specified as Uq = 0.05 (in 
lattice units). 
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Figure 2: Setup for the simulation of flow past square obstacle 



At the inlet (side AC in Fig. EJ), the populations are replaced with the equilib- 
rium values that correspond to the inlet velocity and density. As, the simulation 
result is not very sensitive to the exact condition specified at the inlet this simple 
lowest order approximation at the inlet is sufficient. However, the simulation is very 
sensitive to the condition employed at the exit (side BD in Fig. which should 
be provided in a way that does not affect the bulk flow strongly. The sensitivity of 
the outlet condition is known for this problem EEl We use the simplest possible 
strategy of taking a domain long enough in the direction of flow (up to 45 -L c h ar ) ■ 
At the exit, the populations pointing towards the domain are simply replaced by 
the equilibrium values that correspond to the extrapolated velocity and density. On 
the top and bottom surfaces (side AB and CD in Fig. (2J, the free-slip boundary 
condition is imposed ^ On the cylinder walls (square qrst in Fig. 0), the diffusive 
boundary condition described in the previous section is used. 

The analysis of the unsteady data is performed using the discrete Fourier trans- 
form. The vortex shedding frequency, n, is obtained from the discrete Fourier trans- 
form (in time) of the instantaneous velocity at a monitoring point. The unsteady 
data is collected from two or three point downstream of the square cylinder. 

4. A comparison of Experiment and Simulation 

Numerical simulations of the unsteady flow in the wake of rectangular cylin- 
ders immersed in an infinite stream have been carried out by many authors in the 
past. Much of the effort concentrated either in studying the behavior for moderate 
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Reynolds number flows I I I I I I or verifying sub-gri d models of t urbulence using 
experimental data at high Reynolds numbei|2SI 2a ^ l 31 l a2 l a3 l 34 l a5 l (Re = 21400). 

In this work, we will investigate the behavior of Strouhal number as a function 
of Reynolds number for moderate as well as high Reynolds number flows. Experi- 
mentally, two distinct regimes are known for this flow 1 * 1 ' 1 . First, a sharp transition 
in the value of the Strouhal number is observed for Reynolds number, Re ~ 70. 
Further, for high Reynolds number flows (Re S> 1000), the Strouhal number is al- 
most independence of the Reynolds number and attainds an almost constant value 
of 0.13. We will show that indeed using entropic lattice Boltzmann method the 
experimental behavior can be predicted in a quantitative fashion. 

In Figures |21 and 01 show the streamlines for Re = 100 and Re = 3000, respec- 
tively, after the initial transient dies out. As it can be seen from the plot, at high 
Reynolds number the numerical simulation is not well resolved. In addition, the 
errors from the outflow boundary starts to corrupt the simulation. 

Fig. 0shows a comparison of the experimental ^Hl and simulation results. Table 
H summarizes our simulation results with some of the earlier reported works. It is 
clear that the present method is able to capture the transition around Re ~ 70 as 
well as the asymptotic behavior at high Reynolds numbers in a quantitatively cor- 
rect fashion. It is worth mentioning that the resolution used for the high Reynolds 
number simulations (Re 3> 250) is quite low in comparison to the standard lattice 
Boltzmann method. In order to numerically simulate this flow for Re 3> 1000, using 
a standard lattice Boltzmann method the required number of grid-points would be 
at least O(10 7 ) to avoid numerical instabilities. The fact that an underresolved sim- 
ulation (number of grid-points used is O(10 5 )), without any explicit sub-grid scale 



Entropic lattice Simulation of the flow past square 



700 
600 
500 
400 
300 
200 

200 300 400 500 600 700 800 
X 

Figure 4: A snap-shot of streamfunction at Re = 3000 with L c h a r — 24. 

model, can predict the experimentally observed variation of the Strouhal number 
(a coarse-feature of the flow) indicates that the present method is a candidate for 
efficient coarse-grained simulations of high Reynolds numbers flows. Similar results 
were obtained in the case of homogeneous decaying turbulence 

Finally, at very-high Reynolds number Re = 20000, with a grid of L c h ar = 32 
and L = AbL c har the Strouhal number is fluctuating between 0.09 — 0.15, and errors 
from the boundary corrupted the simulation before the Strouhal number can reach 
a final steady state. 



Table 1. Variation of Strouhal number as a function of Re. 





Re 


Experiment I^EI 


ELBM 


LBM 


fdIU 


F^ 




37 


N. A. 


0.098 


N. A. 


N. A. 


N. A. 




81 


0.115-0.130 


0.122 


0.132 




0.132 




110 


0.143-0.145 


0.141 


0.15 




N. A. 




250 


0.140-0.143 


0.137 




0.16-0.17 


0.154 




500 


0.125-0.127 








0.174 




1000 


0.122-0.123 


0.123 




0.14-0.15 


N. A. 




2800 


0.128-0.129 




N.A 


0.14-0.15 


N. A. 




3000 


0.128-0.129 


0.133 


N.A 


N.A. 


N. A. 




5000 


0.127-0.131 


0.134 


N.A 


N.A. 


N. A. 




20000 


0.13-0.134 


0.09-0.15 


N.A 


N.A. 


N. A. 



Before concluding we will like to comment on the importance of the dimension- 
ality of the flow. Above some critical Reynolds number, the flow becomes three- 
dimensional. Nevertheless, our two-dimensional simulations accurately capture the 
Strouhal number dependence on the Reynolds number. We are currently extending 
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Figure 5: A comparison of the present simulation with the available experimental 
result (Ref . 26) . The figure shows the variation of the Strouhal number as a function 
of Reynolds number. The simulation points are shown as rectangles with a filled 
circle. 
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the code to study the three-dimensional effects in detail. 
5. Conclusions 

We have shown that ELBM can quantitatively predict the variation of vortex shed- 
ding frequency as a function of Reynolds number even for Re » 1000 for the flow 
past a square cylinder. This shows that at least large scale features of the motion 
are effectively captured by the built-in sub-grid scale model of the ELBM. Further 
detailed study using three-dimensional model are needed to check the accuracy and 
efficiency of the method for Re ~ 0(10000). 
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